home *** CD-ROM | disk | FTP | other *** search
/ Amiga Games: Greatest Hits 1996 / Amiga Games: Greatest Hits 1996.iso / archive / userbox / publicdomain / fractal.lha / fractal / source_amiga / enc.c < prev    next >
C/C++ Source or Header  |  1996-07-02  |  40KB  |  960 lines
  1. /**************************************************************************/
  2. /* Encode a byte image using a fractal scheme with a quadtree partition   */
  3. /*                                                                        */
  4. /*       Copyright 1993,1994 Yuval Fisher. All rights reserved.           */
  5. /*                                                                        */
  6. /* Version 0.03 3/14/94                                                   */
  7. /**************************************************************************/
  8.  
  9. /* see CHANGES.AMIGA for which changes done for the Amiga port
  10.    by Andreas R. Kleinert
  11.  
  12.    V1.0: 02-July-96
  13. */
  14.  
  15. #include <stdio.h>
  16. #include <math.h>
  17.  
  18. /* added by ak */
  19. #include <string.h>
  20. #include <stdlib.h>
  21. #include <m68881.h>
  22. /* eo - added by ak */
  23.  
  24. #define DEBUG 0
  25. #define GREY_LEVELS 255
  26.  
  27. #define bound(a)   ((a) < 0.0 ? 0 : ((a)>255.0? 255 : a))
  28. #define IMAGE_TYPE unsigned char /* may be different in some applications */
  29.  
  30. /* various function declarations to keep compiler warnings away. ANSI     */
  31. /* prototypes can go here, for the hearty.                                */
  32. void fatal();
  33. #ifdef EXCLUDED_AK
  34. char *malloc();
  35. char *strcpy(); 
  36. #endif /* EXCLUDED_AK */
  37.  
  38. /* The following #define allocates an hsize x vsize  matrix of type TYPE */
  39. #define matrix_allocate(matrix, hsize, vsize, TYPE) {\
  40.     TYPE *imptr; \
  41.     int _i; \
  42.     matrix = (TYPE **)malloc((vsize)*sizeof(TYPE *));\
  43.     imptr = (TYPE*)malloc((long)(hsize)*(long)(vsize)*sizeof(TYPE));\
  44.     if (imptr == NULL) \
  45.         fatal("\nNo memory in matrix allocate."); \
  46.     for (_i = 0; _i<vsize; ++_i, imptr += hsize) \
  47.         matrix[_i] = imptr; \
  48.  }
  49.  
  50. #define swap(a,b,TYPE)           {TYPE _temp; _temp=b; b=a; a= _temp;}
  51.  
  52. IMAGE_TYPE **image;          /* The input image data                      */ 
  53. double     **domimage[4];    /* Decimated input image used for domains    */
  54.  
  55. double max_scale = 1.0;      /* Maximum allowable grey level scale factor */
  56.  
  57. int    s_bits = 5,           /* Number of bits used to store scale factor */
  58.        o_bits = 7,           /* Number of bits used to store offset       */
  59.        min_part = 4,         /* Min and max _part determine a range of    */
  60.        max_part = 6,         /* Range sizes from hsize>>min to hsize>>max */
  61.        dom_step = 1,         /* Density of domains relative to size       */
  62.        dom_step_type = 0,    /* Flag for dom_step a multiplier or divisor */
  63.        dom_type = 0,         /* Method of generating domain pool 0,1,2..  */
  64.        only_positive = 0,    /* A flag specifying use of positive scaling */
  65.        subclass_search = 0,  /* A flag specifying classes searched        */
  66.        fullclass_search = 0, /* A flag specifying classes searched        */
  67.        *bits_needed,         /* Number of bits to encode domain position. */
  68.        zero_ialpha,          /* The const ialpha when alpha = 0           */
  69.        max_exponent;         /* The max power of 2 side of square image   */
  70.                              /* that fits in our input image.             */
  71.  
  72.  
  73.                              /* The class_transform gives the transforms  */
  74.                              /* between classification numbers for        */
  75.                              /* negative scaling values, when brightest   */
  76.                              /* becomes darkest, etc...                   */
  77. int    class_transform[2][24] = {23,17,21,11,15,9,22,16,19,5,13,3,20,10,18,
  78.                                  4,7,1,14,8,12,2,6,0,
  79.                                  16,22,10,20,8,14,17,23,4,18,2,12,11,21,5,
  80.                                  19,0,6,9,15,3,13,1,7};
  81.  
  82.                              /* rot_transform gives the rotations for     */
  83.                              /* domains with negative scalings.           */
  84. int    rot_transform[2][8] = {7,4,5,6,1,2,3,0, 2,3,0,1,6,7,4,5};
  85.  
  86. struct domain_data {
  87.        int *no_h_domains,    /* The number of domains horizontally for    */
  88.            *no_v_domains,    /* each size.                                */
  89.            *domain_hsize,    /* The size of the domain.                   */
  90.            *domain_vsize,    /* The size of the domain.                   */
  91.            *domain_hstep,    /* The density of the domains.               */
  92.            *domain_vstep;    /* The density of the domains.               */
  93. struct domain_pixels {       /* This is a three (sigh) index array that   */
  94.        int dom_x, dom_y;     /* dynamically allocated. The first index is */
  95.        double sum,sum2;      /* the domain size, the other are two its    */
  96.        int sym;              /* position. It contains the sum and sum^2   */
  97. } ***pixel;                  /* of the pixel values in the domains, which */
  98. } domain;                    /* are computed just once.                   */ 
  99.                              
  100.  
  101. struct classified_domain {             /* This is a list which containes  */
  102.        struct domain_pixels *the;      /* pointers to  the domain data    */
  103.        struct classified_domain *next; /* in the structure above. There   */
  104. } **the_domain[3][24];                 /* are three classes with 24 sub-  */
  105.                                        /* classes. Using this array, only */
  106.                                        /* domains and ranges in the same  */
  107.                                        /* class are compared..            */
  108.                                        /* The first pointer points to the */
  109.                                        /* domain size the the second to   */
  110.                                        /* list of domains.                */
  111.  
  112. FILE *output;                /* Output FILE containing compressed data    */
  113.  
  114. main(argc,argv)
  115. int argc;
  116. char **argv;
  117. /* Usage: quadfrac [tol [inputfilename [outputfilename [hsize [vsize]]]]]  */
  118. {
  119.     /* Defaults are set initially */
  120.     double         tol = 8.0;            /* Tolerance value for quadtree.  */
  121.     char           inputfilename[200];
  122.     char           outputfilename[200];
  123.     int            i,j,k,
  124.                    hsize = -1,           /* The horizontal and vertical    */
  125.                    vsize = -1;           /* size of the input image.       */
  126.     long           stripchar=0;          /* chars to ignore in input file. */
  127.     FILE           *input; 
  128.  
  129.     inputfilename[0] = 1;  /* We initially set the input to this and */
  130.     outputfilename[0] = 1; /* then check if the input/output names   */
  131.                            /* have been set below.                   */
  132.  
  133.     /* scan through the input line and read in the arguments */
  134.     for (i=1; i<argc; ++i) 
  135.        if (argv[i][0] != '-' ) 
  136.             if (inputfilename[0] == 1)
  137.                 strcpy(inputfilename, argv[i]);
  138.             else if (outputfilename[0] == 1)
  139.                 strcpy(outputfilename, argv[i]);
  140.             else;
  141.        else { /* we have a flag */
  142.             if (strlen(argv[i]) == 1) break;
  143.             switch(argv[i][1]) {
  144.                case 't': tol = atof(argv[++i]);
  145.                          break;
  146.                case 'S': stripchar = atoi(argv[++i]);
  147.                          break;
  148.                case 'x': 
  149.                case 'w': hsize = atoi(argv[++i]);
  150.                          break;
  151.                case 'y': 
  152.                case 'h': vsize = atoi(argv[++i]);
  153.                          break;
  154.                case 'D': dom_type = atoi(argv[++i]);
  155.                          break;
  156.                case 'd': dom_step = atoi(argv[++i]);
  157.                          if (dom_step < 0 || dom_step > 15) 
  158.                             fatal("\n Bad domain step.");
  159.                          break;
  160.                case 's': s_bits = atoi(argv[++i]);
  161.                          break;
  162.                case 'o': o_bits = atoi(argv[++i]);
  163.                          break;
  164.                case 'm': min_part = atoi(argv[++i]);
  165.                          break;
  166.                case 'M': max_part = atoi(argv[++i]);
  167.                          break;
  168.                case 'e': dom_step_type= 1;
  169.                          break;
  170.                case 'p': only_positive = 1;
  171.                          break;
  172.                case 'f': subclass_search = 1;
  173.                          break;
  174.                case 'F': fullclass_search = 1;
  175.                          break;
  176.                case 'N': max_scale = atof(argv[++i]);
  177.                          break;
  178.                case '?':
  179.                case 'H':
  180.                default:
  181.            printf("\nUsage: enc -[options] [inputfile [outputfile]]");
  182.            printf("\nOptions are: (# = number), defaults show in ()");
  183.            printf("\n  -t # tolerance criterion for fidelity. (%lf)", tol);
  184.            printf("\n  -m # minimum quadtree partitions.  (%d)",min_part);
  185.            printf("\n  -M # maximum quadtree partitions. (%d)",max_part);
  186.            printf("\n  -S # number of input bytes to ignore. (%ld)",stripchar);
  187.            printf("\n  -w # width (horizontal size) of input data. (256)");
  188.            printf("\n  -h # height (vertical size) of input data. (256)");
  189.            printf("\n  -d # domain step size. (%d)", dom_step);
  190.            printf("\n  -D # method {0,1,2} for domain pool (%d)",dom_type);
  191.            printf("\n  -s # number of scaling quantizing bits. (%d)",s_bits);
  192.            printf("\n  -o # number of offset quantizing bits. (%d)",o_bits);
  193.            printf("\n  -N # maximum scaling in encoding. (%lf)",max_scale);
  194.            printf("\n  -e   domain step as multiplier not divisor. (off)");
  195.            printf("\n  -p   use only positive scaling (for speed). (off)");
  196.            printf("\n  -f   search 24 domain classes (for fidelity). (off)");
  197.            printf("\n  -F   search 3 domain classes (for fidelity). (off)");
  198.            fatal("\n       -F and -f can be used together.");
  199.             }
  200.        }
  201.    
  202.     if (inputfilename[0] == 1) strcpy(inputfilename, "lenna.dat");
  203.     if (outputfilename[0] == 1) strcpy(outputfilename, "lenna.trn");
  204.  
  205.     if (hsize == -1) 
  206.         if (vsize == -1) hsize = vsize = 256;
  207.         else hsize = vsize;
  208.     else 
  209.         if (vsize == -1) vsize = hsize;
  210.  
  211.     /* allocate memory for the input image. Allocating one chunck saves  */
  212.     /* work and time later.                                              */
  213.     matrix_allocate(image, hsize, vsize, IMAGE_TYPE)
  214.     matrix_allocate(domimage[0], hsize/2, vsize/2, double)
  215.     matrix_allocate(domimage[1], hsize/2, vsize/2, double)
  216.     matrix_allocate(domimage[2], hsize/2, vsize/2, double)
  217.     matrix_allocate(domimage[3], hsize/2, vsize/2, double)
  218.  
  219.     /* max_ & min_ part are variable, so this must be run time allocated */
  220.     bits_needed = (int *)malloc(sizeof(int)*(1+max_part-min_part));
  221.  
  222.     if ((input = fopen(inputfilename, "r")) == NULL)
  223.         fatal("Can't open input file.");
  224.  
  225.     /* skip the first  stripchar  chars */ 
  226.     fseek(input, stripchar, 0); 
  227.     i = fread(image[0], sizeof(IMAGE_TYPE), hsize*vsize, input);
  228.     fclose(input);
  229.  
  230.     if (i < hsize*vsize) 
  231.         fatal("Not enough input data in the input file.");
  232.     else
  233.         printf("%dx%d=%d pixels read from %s.", hsize,vsize,i,inputfilename);
  234.  
  235.     /* allcate memory for domain data and initialize it */
  236.     compute_sums(hsize,vsize);
  237.  
  238.     if ((output = fopen(outputfilename, "w")) == NULL)
  239.         fatal("Can't open output file.");
  240.  
  241.     /* output some data into the outputfile.                       */
  242.     pack(4,(long)min_part,output);
  243.     pack(4,(long)max_part,output);
  244.     pack(4,(long)dom_step,output);
  245.     pack(1,(long)dom_step_type,output);
  246.     pack(2,(long)dom_type,output);
  247.     pack(12,(long)hsize,output);
  248.     pack(12,(long)vsize,output);
  249.  
  250.     /* This is the quantized value of zero scaling.. needed later */
  251.     zero_ialpha = 0.5 + (max_scale)/(2.0*max_scale)*(1<<s_bits);
  252.  
  253.     /* The following routine takes a rectangular image and calls the */
  254.     /* quadtree routine to encode square sum-images in it.           */
  255.     /* the tolerance is a parameter since in some applications different */
  256.     /* regions of the image may need to be compressed to different tol's */
  257.     printf("\nEncoding Image.....");
  258.     fflush(stdout);
  259.     partition_image(0, 0, hsize,vsize, tol);
  260.     printf("Done.");
  261.     fflush(stdout);
  262.  
  263.     /* stuff the last byte if needed */
  264.     pack(-1,(long)0,output);
  265.     
  266.     fclose(output);
  267.     i = pack(-2,(long)0,output);
  268.     printf("\n Compression = %lf from %d bytes written in %s.\n",
  269.             (double)(hsize*vsize)/(double)i, i, outputfilename);
  270.  
  271.     /* Free allocated memory*/
  272.     free(bits_needed);
  273.     free(domimage[0]);
  274.     free(domimage[1]);
  275.     free(domimage[2]);
  276.     free(domimage[3]);
  277.     free(domain.no_h_domains);
  278.     free(domain.no_v_domains);
  279.     free(domain.domain_hsize);
  280.     free(domain.domain_vsize);
  281.     free(domain.domain_hstep);
  282.     free(domain.domain_vstep);
  283.     for (i=0; i <= max_part-min_part; ++i) 
  284.         free(domain.pixel[i]);
  285.     free(domain.pixel);
  286.     free(image[0]);
  287.     for (i=0; i <= max_part-min_part; ++i) 
  288.     for (k=0; k<3; ++k) 
  289.      for (j=0; j<24; ++j) list_free(the_domain[k][j][i]);
  290.     return(0);
  291. }
  292.  
  293. /* ************************************************************** */
  294. /* free memory allocated in the list structure the_domain         */
  295. /* ************************************************************** */
  296. list_free(node)
  297. struct classified_domain *node;
  298. {
  299.      if (node->next != NULL)
  300.        list_free(node->next);
  301.      free(node);
  302. }
  303.  
  304. /* ************************************************************** */
  305. /* return the average pixel value of a region of the image.       */
  306. /* ************************************************************** */
  307. void average(x,y,xsize,ysize, psum, psum2)
  308. int x,y,xsize,ysize;
  309. double *psum, *psum2;
  310. {
  311.     register int i,j,k;
  312.     register double pixel;
  313.     *psum = *psum2 = 0.0;
  314.     k = ((x%2)<<1) + y%2;
  315.     x >>= 1; y >>= 1;
  316.     xsize >>= 1; ysize >>= 1;
  317.     for (i=x; i<x+xsize; ++i)
  318.     for (j=y; j<y+ysize; ++j) {
  319.        pixel = domimage[k][j][i];
  320.        *psum += pixel;
  321.        *psum2 += pixel*pixel;
  322.     }
  323. }
  324.  
  325. /* ************************************************************** */
  326. /* return the average pixel value of a region of the image. This  */
  327. /* routine differs from the previous in one slight way. It does   */
  328. /* not average 2x2 sub-images to pixels. This is needed for clas- */
  329. /* sifying ranges rather than domain where decimation is needed.  */
  330. /* ************************************************************** */
  331. void average1(x,y,xsize,ysize, psum, psum2)
  332. int x,y,xsize,ysize;
  333. double *psum, *psum2;
  334. {
  335.     register int i,j;
  336.     register double pixel;
  337.     *psum = *psum2 = 0.0;
  338.  
  339.     for (i=x; i<x+xsize; ++i)
  340.     for (j=y; j<y+ysize; ++j) {
  341.        pixel = (double)image[j][i];
  342.        *psum += pixel;
  343.        *psum2 += pixel*pixel;
  344.     }
  345. }
  346.  
  347. /* ************************************************************** */
  348. /* Take a region of the image at x,y and classify it.             */
  349. /* The four quadrants of the region are ordered from brightest to */
  350. /* least bright average value, then it is rotated into one of the */
  351. /* three cannonical orientations possible with the brightest quad */
  352. /* in the upper left corner.                                      */
  353. /* The routine returns two indices that are class numbers: pfirst */
  354. /* and psecond; the symmetry operation that bring the square into */
  355. /* cannonical position; and the sum and sum^2 of the pixel values */
  356. /* ************************************************************** */
  357. classify(x, y, xsize, ysize, pfirst, psecond, psym, psum, psum2, type)
  358. int x,y,xsize,ysize,   /* position, size of subimage to be classified */
  359.     *pfirst, *psecond, /* returned first and second class numbers     */
  360.     *psym;             /* returned symmetry operation that brings the */
  361.                        /* subimage to cannonical position.            */
  362. double *psum, *psum2;  /* returned sum and sum^2 of pixel values      */
  363. int type;              /* flag for decimating (for domains) or not    */ 
  364. {
  365.  
  366.     int order[4], i,j; 
  367.     double a[4],a2[4];
  368.     void (*average_func)();
  369.     
  370.     if (type == 2) average_func = average; else average_func = average1;
  371.  
  372.     /* get the average values of each quadrant                         */
  373.  
  374.  
  375.     (*average_func)(x,y,                 xsize/2,ysize/2,   &a[0], &a2[0]);
  376.     (*average_func)(x,y+ysize/2,         xsize/2,ysize/2,   &a[1], &a2[1]);
  377.     (*average_func)(x+xsize/2,y+ysize/2, xsize/2,ysize/2,   &a[2], &a2[2]);
  378.     (*average_func)(x+xsize/2,y,         xsize/2,ysize/2,   &a[3], &a2[3]);
  379.  
  380.     *psum = a[0] + a[1] + a[2] + a[3];
  381.     *psum2 =  a2[0] + a2[1] + a2[2] + a2[3];
  382.  
  383.     for (i=0; i<4; ++i) {
  384.          /* after the sorting below order[i] is the i-th brightest   */
  385.          /* quadrant.                                                */
  386.          order[i] = i; 
  387.          /* convert a2[] to store the variance of each quadrant      */
  388.          a2[i] -= (double)(1<<(2*type))*a[i]*a[i]/(double)(xsize*ysize);
  389.     }
  390.  
  391.     /* Now order the average value and also in order[],   which will */
  392.     /* then tell us the indices (in a[]) of the brightest to darkest */
  393.     for (i=2; i>=0; --i)
  394.     for (j=0; j<=i; ++j)
  395.         if (a[j]<a[j+1]) {
  396.             swap(order[j], order[j+1],int)
  397.             swap(a[j], a[j+1],double)
  398.     }
  399.  
  400.     /* because of the way we ordered the a[] the rotation can be */
  401.     /* read right off of order[]. That will make the brightest   */
  402.     /* quadrant be in the upper left corner. But we must still   */ 
  403.     /* decide which cannonical class the image portion belogs    */
  404.     /* to and whether to do a flip or just a rotation. This is   */
  405.     /* the following table summarizes the horrid lines below     */
  406.     /* order      class            do a rotation                 */
  407.     /* 0,2,1,3      0                   0                        */
  408.     /* 0,2,3,1      0                   1                        */
  409.     /* 0,1,2,3      1                   0                        */
  410.     /* 0,3,2,1      1                   1                        */
  411.     /* 0,1,3,2      2                   0                        */
  412.     /* 0,3,1,2      2                   1                        */
  413.  
  414.     *psym = order[0]; 
  415.     /* rotate the values */
  416.     for (i=0; i<4; ++i)
  417.         order[i] = (order[i] - (*psym) + 4)%4; 
  418.  
  419.     for (i=0; order[i] != 2; ++i); 
  420.     *pfirst = i-1;
  421.     if (order[3] == 1 || (*pfirst == 2 && order[2] == 1)) *psym += 4;
  422.  
  423.     /* Now further classify the sub-image by the variance of its    */
  424.     /* quadrants. This give 24 subclasses for each of the 3 classes */
  425.     for (i=0; i<4; ++i) order[i] = i; 
  426.  
  427.     for (i=2; i>=0; --i)
  428.     for (j=0; j<=i; ++j)
  429.         if (a2[j]<a2[j+1]) {
  430.             swap(order[j], order[j+1],int)
  431.             swap(a2[j], a2[j+1],double)
  432.     }
  433.  
  434.     /* Now do the symmetry operation */
  435.     for (i=0; i<4; ++i)
  436.         order[i] = (order[i] - (*psym%4) + 4)%4; 
  437.     if (*psym > 3) 
  438.         for (i=0; i<4; ++i)
  439.            if (order[i]%2) order[i] = (2 + order[i])%4;
  440.  
  441.     /* We want to return a class number from 0 to 23 depending on */
  442.     /* the ordering of the quadrants according to their variance  */
  443.     *psecond = 0;
  444.     for (i=2; i>=0; --i)
  445.     for (j=0; j<=i; ++j)
  446.         if (order[j] > order[j+1]) {
  447.             swap(order[j],order[j+1], int);
  448.             if (order[j] == 0 || order [j+1] == 0) 
  449.                  *psecond += 6;
  450.             else if (order[j] == 1 || order [j+1] == 1) 
  451.                  *psecond += 2;
  452.             else if (order[j] == 2 || order [j+1] == 2) 
  453.                  *psecond += 1;
  454.         }
  455. }
  456.  
  457. /* ************************************************************ */
  458. /* Compute sum and sum^2 of pixel values in domains for use in  */
  459. /* the rms computation later. Since a domain is compared with   */
  460. /* many ranges, doing this just once saves a lot of computation */
  461. /* This routine also fills a list structure with the domains    */
  462. /* as they are classified and creates the memory for the domain */
  463. /* data in a matrix.                                            */
  464. /* ************************************************************ */
  465. compute_sums(hsize,vsize)
  466. int hsize,vsize;
  467. {
  468.     int i,j,k,l,
  469.          domain_x, 
  470.          domain_y, 
  471.          first_class, 
  472.          second_class,
  473.          domain_size, 
  474.          domain_step_size, 
  475.          size,  
  476.          x_exponent, 
  477.          y_exponent;
  478.  
  479.     struct classified_domain *node;
  480.  
  481.     printf("\nComputing domain sums... ");
  482.     fflush(stdout); 
  483.  
  484.     /* pre-decimate the image into domimage to avoid having to      */
  485.     /* do repeated averaging of 2x2 pixel groups.                   */
  486.     /* There are 4 ways to decimate the image, depending on the     */
  487.     /* location of the domain, odd or even address.                 */
  488.     for (i=0; i<2; ++i)
  489.     for (j=0; j<2; ++j)
  490.     for (k=i; k<hsize-i; k += 2)
  491.     for (l=j; l<vsize-j; l += 2) 
  492.         domimage[(i<<1)+j][l>>1][k>>1] = 
  493.                      ((double)image[l][k] + (double)image[l+1][k+1] +
  494.                      (double)image[l][k+1] + (double)image[l+1][k])*0.25;
  495.  
  496.  
  497.     /* Allocate memory for the sum and sum^2 of domain pixels        */
  498.     /* We first compute the size of the largest square that fits in  */
  499.     /* the image.                                                    */
  500.     x_exponent = (int)floor(log((double)hsize)/log(2.0));
  501.     y_exponent = (int)floor(log((double)vsize)/log(2.0));
  502.    
  503.     /* exponent is min of x_ and y_ exponent */
  504.     max_exponent = (x_exponent > y_exponent ? y_exponent : x_exponent);
  505.  
  506.     /* size is the size of the largest square that fits in the image */
  507.     /* It is used to compute the domain and range sizes.             */
  508.     size =  1<<max_exponent; 
  509.  
  510.     if (max_exponent < max_part)
  511.       fatal("Reduce maximum number of quadtree partitions.");
  512.     if (max_exponent-2 < max_part)
  513.       printf("\nWarning: so many quadtree partitions yield absurd ranges.");
  514.  
  515.     i = max_part - min_part + 1;
  516.     domain.no_h_domains = (int *)malloc(sizeof(int)*i);
  517.     domain.no_v_domains = (int *)malloc(sizeof(int)*i);
  518.     domain.domain_hsize = (int *)malloc(sizeof(int)*i);
  519.     domain.domain_vsize = (int *)malloc(sizeof(int)*i);
  520.     domain.domain_hstep = (int *)malloc(sizeof(int)*i);
  521.     domain.domain_vstep = (int *)malloc(sizeof(int)*i);
  522.  
  523.     domain.pixel= (struct domain_pixels ***)
  524.               malloc(i*sizeof(struct domain_pixels **));
  525.     if (domain.pixel == NULL) fatal("No memory for domain pixel sums.");
  526.  
  527.     for (i=0; i <= max_part-min_part; ++i) {
  528.        /* first compute how many domains there are horizontally */
  529.        domain.domain_hsize[i] = size >> (min_part+i-1);
  530.        if (dom_type == 2) 
  531.              domain.domain_hstep[i] = dom_step; 
  532.        else if (dom_type == 1)
  533.              if (dom_step_type == 1) 
  534.                domain.domain_hstep[i] = (size >> (max_part - i-1))*dom_step;
  535.              else 
  536.                domain.domain_hstep[i] = (size >> (max_part - i-1))/dom_step;
  537.        else 
  538.              if (dom_step_type == 1) 
  539.                domain.domain_hstep[i] = domain.domain_hsize[i]*dom_step;
  540.              else 
  541.                domain.domain_hstep[i] = domain.domain_hsize[i]/dom_step;
  542.  
  543.        domain.no_h_domains[i] = 1+(hsize-domain.domain_hsize[i])/
  544.                                                   domain.domain_hstep[i];
  545.        /* bits_needed[i][0] = ceil(log(domain.no_h_domains[i])/log(2.0));  */
  546.  
  547.        /* now compute how many domains there are vertically. The sizes  */
  548.        /* are the same for square domains, but not for rectangular ones */
  549.        domain.domain_vsize[i] = size >> (min_part+i-1);
  550.        if (dom_type == 2) 
  551.              domain.domain_vstep[i] = dom_step; 
  552.        else if (dom_type == 1)
  553.              if (dom_step_type == 1) 
  554.                domain.domain_vstep[i] = (size >> (max_part - i-1))*dom_step;
  555.              else
  556.                domain.domain_vstep[i] = (size >> (max_part - i-1))/dom_step;
  557.        else 
  558.              if (dom_step_type == 1) 
  559.                domain.domain_vstep[i] = domain.domain_vsize[i]*dom_step;
  560.              else
  561.                domain.domain_vstep[i] = domain.domain_vsize[i]/dom_step;
  562.  
  563.        domain.no_v_domains[i] = 1+(vsize-domain.domain_vsize[i])/
  564.                                                   domain.domain_vstep[i];
  565.  
  566.        /* Now compute the number of bits needed to store the domain data */
  567.        bits_needed[i] = ceil(log((double)domain.no_h_domains[i]*
  568.                                  (double)domain.no_v_domains[i])/log(2.0)); 
  569.  
  570.        matrix_allocate(domain.pixel[i], domain.no_h_domains[i],
  571.                      domain.no_v_domains[i], struct domain_pixels) 
  572.  
  573.     }
  574.  
  575.     /* allocate and zero the list containing the classified domain data */
  576.     i = max_part - min_part + 1;
  577.     for (first_class = 0; first_class < 3; ++first_class)
  578.     for (second_class = 0; second_class < 24; ++second_class) {
  579.        the_domain[first_class][second_class] = 
  580.                            (struct classified_domain **) 
  581.                            malloc(i*sizeof(struct classified_domain *));
  582.        for (j=0; j<i; ++j)
  583.                     the_domain[first_class][second_class][j] = NULL;
  584.     }
  585.  
  586.     /* Precompute sum and sum of squares for domains                 */
  587.     /* This part can get faster for overlapping domains if repeated  */
  588.     /* sums are avoided                                              */
  589.     for (i=0; i <= max_part-min_part; ++i) {
  590.       for (j=0,domain_x=0; j<domain.no_h_domains[i]; ++j,
  591.              domain_x+=domain.domain_hstep[i]) 
  592.       for (k=0,domain_y=0; k<domain.no_v_domains[i]; ++k,
  593.              domain_y+=domain.domain_vstep[i]) {
  594.         classify(domain_x, domain_y, 
  595.                  domain.domain_hsize[i], 
  596.                  domain.domain_vsize[i], 
  597.                      &first_class, &second_class,
  598.                      &domain.pixel[i][k][j].sym, 
  599.                      &domain.pixel[i][k][j].sum, 
  600.                      &domain.pixel[i][k][j].sum2, 2);
  601.  
  602.         /* When the domain data is referenced from the list, we need to */
  603.         /* know where the domain is.. so we have to store the position  */
  604.         domain.pixel[i][k][j].dom_x = j;
  605.         domain.pixel[i][k][j].dom_y = k;
  606.         node = (struct classified_domain *)
  607.                                 malloc(sizeof(struct classified_domain));
  608.         
  609.         /* put this domain in the classified list structure */
  610.         node->the = &domain.pixel[i][k][j];
  611.         node->next = the_domain[first_class][second_class][i];
  612.         the_domain[first_class][second_class][i] = node;
  613.       }
  614.     }
  615.  
  616.     /* Now we make sure no domain class is actually empty.  */
  617.     for (i=0; i <= max_part-min_part; ++i) 
  618.     for (first_class = 0; first_class < 3; ++first_class)
  619.     for (second_class = 0; second_class < 24; ++second_class) 
  620.        if (the_domain[first_class][second_class][i] == NULL) {
  621.            node = (struct classified_domain *)
  622.                                 malloc(sizeof(struct classified_domain));
  623.            node->the = &domain.pixel[i][0][0];
  624.            node->next = NULL;
  625.            the_domain[first_class][second_class][i] = node;
  626.        }
  627.  
  628.     printf("Done.");
  629.     fflush(stdout); 
  630. }
  631.  
  632. /* ************************************************************ */
  633. /* pack value using size bits and output into foutf */
  634. /* ************************************************************ */
  635. int pack(size, value, foutf)
  636. int size; long int value;
  637. FILE *foutf;
  638. {
  639.      int i;
  640.      static int ptr = 1, /* how many bits are packed in sum so far */
  641.                 sum = 0, /* packed bits */ 
  642.                 num_of_packed_bytes = 0; /* total bytes written out */
  643.  
  644.     /* size == -1 means we are at the end, so write out what is left */
  645.     if (size == -1 && ptr != 1) {
  646.         fputc(sum<<(8-ptr), foutf);
  647.         ++num_of_packed_bytes;
  648.         return(0);
  649.     }
  650.     
  651.     /* size == -2 means we want to know how many bytes we have written */
  652.     if (size == -2) 
  653.             return(num_of_packed_bytes);
  654.     
  655.     for (i=0; i<size; ++i, ++ptr, value = value>>1, sum = sum<<1) {
  656.         if (value & 1) sum |= 1;
  657.  
  658.          if (ptr == 8) {
  659.             fputc(sum,foutf);
  660.             ++num_of_packed_bytes;
  661.             sum=0;
  662.             ptr=0;
  663.         }
  664.      }
  665. }
  666.  
  667. /* ************************************************************ */
  668. /* Compare a range to a domain and return rms and the quantized */
  669. /* scaling and offset values (pialhpa, pibeta).                 */
  670. /* ************************************************************ */
  671. double compare(atx,aty, xsize, ysize, depth, rsum, rsum2, dom_x,dom_y, 
  672.                                   sym_op, pialpha,pibeta)
  673. int atx, aty, xsize, ysize, depth, dom_x, dom_y, sym_op, *pialpha, *pibeta;
  674. double rsum, rsum2;
  675. {
  676.     int i, j, i1, j1, k,
  677.         domain_x, 
  678.         domain_y;        /* The domain position                */
  679.  
  680.     double pixel,
  681.            det,          /* determinant of solution */
  682.            dsum,         /* sum of domain values */
  683.            rdsum = 0,    /* sum of range*domain values   */
  684.            dsum2,        /* sum of domain^2 values   */
  685.            w2 = 0,       /* total number of values tested */
  686.            rms = 0,      /* root means square difference */
  687.            alpha,        /* the scale factor */
  688.            beta;         /* The offset */
  689.  
  690.  
  691.     
  692.     /* xsize = hsize >> depth; */
  693.     /* ysize = vsize >> depth; */
  694.     w2 = xsize * ysize;
  695.  
  696.     dsum = domain.pixel[depth-min_part][dom_y][dom_x].sum;
  697.     dsum2 = domain.pixel[depth-min_part][dom_y][dom_x].sum2;
  698.     domain_x = (dom_x * domain.domain_hstep[depth-min_part]);
  699.     domain_y = (dom_y * domain.domain_vstep[depth-min_part]);
  700.     k = ((domain_x%2)<<1) + domain_y%2;
  701.     domain_x >>= 1;
  702.     domain_y >>= 1;
  703.  
  704.     /* The main statement in this routine is a switch statement which */
  705.     /* scans through the domain and range to compare them. The loop   */
  706.     /* center is the same so we #define it for easy modification      */
  707. #define COMPUTE_LOOP              {                                 \
  708.         pixel = domimage[k][j1][i1];                                \
  709.         rdsum += image[j][i]*pixel;                                 \
  710.         }
  711.     
  712.      switch(sym_op) {
  713.          case 0: for (i=atx, i1 = domain_x; i<atx+xsize; ++i, ++i1)
  714.                  for (j=aty, j1 = domain_y; j<aty+ysize; ++j, ++j1) 
  715.                      COMPUTE_LOOP
  716.                  break;
  717.          case 1: for (j=aty, i1 = domain_x; j<aty+ysize; ++j, ++i1)
  718.                  for (i=atx+xsize-1, j1 = domain_y; i>=atx; --i, ++j1) 
  719.                      COMPUTE_LOOP
  720.                  break;
  721.          case 2: for (i=atx+xsize-1, i1 = domain_x; i>=atx; --i, ++i1)
  722.                  for (j=aty+ysize-1, j1 = domain_y; j>=aty; --j, ++j1) 
  723.                      COMPUTE_LOOP
  724.                  break;
  725.          case 3: for (j=aty+ysize-1, i1 = domain_x; j>=aty; --j, ++i1)
  726.                  for (i=atx, j1 = domain_y; i<atx+xsize; ++i, ++j1) 
  727.                      COMPUTE_LOOP
  728.                  break;
  729.          case 4: for (j=aty, i1 = domain_x; j<aty+ysize; ++j, ++i1)
  730.                  for (i=atx, j1 = domain_y; i<atx+xsize; ++i, ++j1) 
  731.                      COMPUTE_LOOP
  732.                  break;
  733.          case 5: for (i=atx, i1 = domain_x; i<atx+xsize; ++i, ++i1)
  734.                  for (j=aty+ysize-1, j1 = domain_y; j>=aty; --j, ++j1) 
  735.                      COMPUTE_LOOP
  736.                  break;
  737.          case 6: for (j=aty+ysize-1, i1 = domain_x; j>=aty; --j, ++i1)
  738.                  for (i=atx+xsize-1, j1 = domain_y; i>=atx; --i, ++j1) 
  739.                      COMPUTE_LOOP
  740.                  break;
  741.          case 7: for (i=atx+xsize-1, i1 = domain_x; i>=atx; --i, ++i1)
  742.                  for (j=aty, j1 = domain_y; j<aty+ysize; ++j, ++j1) 
  743.                      COMPUTE_LOOP
  744.                  break;
  745.      }
  746.     
  747.      det = (xsize*ysize)*dsum2 - dsum*dsum;
  748.     
  749.      if (det == 0.0) 
  750.          alpha = 0.0; 
  751.       else 
  752.          alpha = (w2*rdsum - rsum*dsum)/det;
  753.         
  754.      if  (only_positive && alpha < 0.0) alpha = 0.0;
  755.      *pialpha = 0.5 + (alpha + max_scale)/(2.0*max_scale)*(1<<s_bits);
  756.      if (*pialpha < 0) *pialpha = 0;
  757.      if (*pialpha >= 1<<s_bits) *pialpha = (1<<s_bits)-1;
  758.  
  759.      /* Now recompute alpha back */
  760.      alpha = (double)*pialpha/(double)(1<<s_bits)*(2.0*max_scale)-max_scale;
  761.         
  762.      /* compute the offset */
  763.      beta = (rsum - alpha*dsum)/w2;
  764.  
  765.      /* Convert beta to an integer */
  766.      /* we use the sign information of alpha to pack efficiently */
  767.      if (alpha > 0.0)  beta += alpha*GREY_LEVELS;
  768.      *pibeta = 0.5 + beta/
  769.              ((1.0+fabs(alpha))*GREY_LEVELS)*((1<<o_bits)-1);
  770.      if (*pibeta< 0) *pibeta = 0;
  771.      if (*pibeta>= 1<<o_bits) *pibeta = (1<<o_bits)-1;
  772.  
  773.      /* Recompute beta from the integer */
  774.      beta = (double)*pibeta/(double)((1<<o_bits)-1)*
  775.                ((1.0+fabs(alpha))*GREY_LEVELS);
  776.      if (alpha > 0.0) beta  -= alpha*GREY_LEVELS;
  777.  
  778.      /* Compute the rms based on the quantized alpha and beta! */
  779.      rms= sqrt((rsum2 + alpha*(alpha*dsum2 - 2.0*rdsum + 2.0*beta*dsum) +
  780.                  beta*(beta*w2 - 2.0*rsum))/w2);
  781.     
  782.      return(rms);
  783. }
  784.  
  785. /* ************************************************************ */
  786. /* Recursively partition an image, computing the best transfoms */
  787. /* ************************************************************ */
  788. quadtree(atx,aty,xsize,ysize,tol,depth)
  789. int atx, aty, xsize, ysize, depth; 
  790. double tol;  /* the tolerance fit  */
  791. {
  792.     int i,
  793.         sym_op,                   /* A symmetry operation of the square */
  794.         ialpha,                   /* Intgerized scalling factor         */
  795.         ibeta,                    /* Intgerized offset                  */
  796.         best_ialpha,              /* best ialpha found                  */
  797.         best_ibeta,
  798.         best_sym_op,
  799.         best_domain_x,
  800.         best_domain_y, 
  801.         first_class, 
  802.         the_first_class, 
  803.         first_class_start,        /* loop beginning and ending values   */
  804.         first_class_end, 
  805.         second_class[2],
  806.         the_second_class, 
  807.         second_class_start,       /* loop beginning and ending values   */
  808.         second_class_end, 
  809.         range_sym_op[2],          /* the operations to bring square to  */
  810.         domain_sym_op;            /* cannonical position.               */
  811.  
  812.     struct classified_domain *node;  /* var for domains we scan through */
  813.  
  814.     double rms, best_rms,          /* rms value and min rms found so far */
  815.            sum=0, sum2=0;          /* sum and sum^2 of range pixels      */
  816.  
  817.  
  818.     /* keep breaking it down until we are small enough */
  819.     if (depth < min_part) {
  820.          quadtree(atx,aty, xsize/2, ysize/2, tol,depth+1);
  821.          quadtree(atx+xsize/2,aty, xsize/2, ysize/2,tol,depth+1);
  822.          quadtree(atx,aty+ysize/2, xsize/2, ysize/2,tol,depth+1);
  823.          quadtree(atx+xsize/2,aty+ysize/2, xsize/2, ysize/2,tol,depth+1);
  824.          return;
  825.     }
  826.  
  827.     /* now search for the best domain-range match and write it out */
  828.     best_rms = 10000000000.0;      /* just a big number */
  829.  
  830.     classify(atx, aty, xsize,ysize, 
  831.               &the_first_class, &the_second_class, 
  832.               &range_sym_op[0], &sum, &sum2, 1);
  833.  
  834.  
  835.     /* sort out how much to search based on -f and -F input flags */
  836.     if (fullclass_search) {
  837.          first_class_start = 0; 
  838.          first_class_end = 3;
  839.     } else {
  840.          first_class_start = the_first_class; 
  841.          first_class_end = the_first_class+1;
  842.     } 
  843.  
  844.     if (subclass_search) {
  845.          second_class_start = 0; 
  846.          second_class_end = 24;
  847.     } else {
  848.          second_class_start = the_second_class; 
  849.          second_class_end = the_second_class+1;
  850.     } 
  851.  
  852.     /* these for loops vary according to the optimization flags we set */
  853.     /* for subclass_search and fullclass_search==1, we search all the  */
  854.     /* domains (except not all rotations).                             */
  855.     for (first_class = first_class_start; 
  856.          first_class < first_class_end; ++first_class)
  857.     for (second_class[0] = second_class_start; 
  858.          second_class[0] < second_class_end; ++second_class[0]) {
  859.  
  860.        /* We must check each domain twice. Once for positive scaling,  */
  861.        /* once for negative scaling. Each has its own class and sym_op */
  862.        if (!only_positive) {
  863.           second_class[1] = 
  864.              class_transform[(first_class == 2 ? 1 : 0)][second_class[0]];
  865.           range_sym_op[1] = 
  866.              rot_transform[(the_first_class == 2 ? 1 : 0)][range_sym_op[0]];
  867.        } 
  868.  
  869.        /* only_positive is 0 or 1, so we may or not scan                */
  870.        for (i=0; i<(2-only_positive); ++i) 
  871.        for (node = the_domain[first_class][second_class[i]][depth-min_part]; 
  872.            node != NULL; 
  873.            node = node->next) {
  874.            domain_sym_op = node->the->sym;
  875.            /* The following if statement figures out how to transform  */
  876.            /* the domain onto the range, given that we know how to get */
  877.            /* each into cannonical position.                           */
  878.            if (((domain_sym_op>3 ? 4: 0) + (range_sym_op[i]>3 ? 4: 0))%8 == 0)
  879.              sym_op = (4 + domain_sym_op%4 - range_sym_op[i]%4)%4;
  880.            else 
  881.              sym_op = (4 + (domain_sym_op%4 + 3*(4-range_sym_op[i]%4))%4)%8;
  882.             
  883.            rms = compare(atx,aty, xsize, ysize,  depth, sum,sum2, 
  884.                                   node->the->dom_x, 
  885.                                   node->the->dom_y, 
  886.                                   sym_op, &ialpha,&ibeta); 
  887.     
  888.            if (rms < best_rms) {
  889.                    best_ialpha = ialpha;
  890.                    best_ibeta = ibeta;
  891.                    best_rms = rms;
  892.                    best_sym_op = sym_op;
  893.                    best_domain_x = node->the->dom_x;
  894.                    best_domain_y = node->the->dom_y;
  895.            }
  896.        }
  897.     }
  898.         
  899.     if (best_rms > tol && depth < max_part) {
  900.        /* We didn't find a good enough fit so quadtree down */
  901.        pack(1,(long)1,output);  /* This bit means we quadtree'd down */
  902.        quadtree(atx,aty, xsize/2, ysize/2, tol,depth+1);
  903.        quadtree(atx+xsize/2,aty, xsize/2, ysize/2, tol,depth+1);
  904.        quadtree(atx,aty+ysize/2, xsize/2, ysize/2, tol,depth+1);
  905.        quadtree(atx+xsize/2,aty+ysize/2, xsize/2, ysize/2, tol,depth+1);
  906.     } else {
  907.        /* The fit was good enough or we just can't get smaller ranges */
  908.        /* So write out the data */
  909.        if (depth < max_part)          /* if we are not at the smallest range */
  910.                pack(1,(long)0,output);/* then we must tell the decoder we    */
  911.                                       /* stopped quadtreeing                 */
  912.        pack(s_bits, (long)best_ialpha, output);
  913.        pack(o_bits, (long)best_ibeta, output);
  914.        /* When the scaling is zero, there is no need to store the domain */
  915.        if (best_ialpha != zero_ialpha) {
  916.           pack(3, (long)best_sym_op, output);
  917.           pack(bits_needed[depth-min_part], (long)(best_domain_y*
  918.             domain.no_h_domains[depth-min_part]+best_domain_x), output);
  919.       }
  920.      }
  921. }
  922.  
  923. /* ************************************************************* */
  924. /* Recursively partition an image, finding the largest contained */
  925. /* square and call the quadtree routine the encode that square.  */
  926. /* This enables us to encode rectangular image easily.           */
  927. /* ************************************************************* */
  928. partition_image(atx, aty, hsize,vsize, tol)
  929. int atx, aty, hsize,vsize;
  930. double tol;
  931. {
  932.    int x_exponent,    /* the largest power of 2 image size that fits */
  933.        y_exponent,    /* horizontally or vertically the rectangle.   */
  934.        exponent,      /* The actual size of image that's encoded.    */
  935.        size, 
  936.        depth; 
  937.    
  938.    x_exponent = (int)floor(log((double)hsize)/log(2.0));
  939.    y_exponent = (int)floor(log((double)vsize)/log(2.0));
  940.    
  941.    /* exponent is min of x_ and y_ exponent */
  942.    exponent = (x_exponent > y_exponent ? y_exponent : x_exponent);
  943.    size = 1<<exponent; 
  944.    depth = max_exponent - exponent;
  945.    quadtree(atx,aty,size,size,tol,depth);
  946.    if (size != hsize) 
  947.       partition_image(atx+size, aty, hsize-size,vsize, tol);
  948.  
  949.    if (size != vsize) 
  950.       partition_image(atx, aty+size, size,vsize-size, tol);
  951. }        
  952.   
  953. /* fatal is for when there is a fatal error... print a message and exit */
  954. void fatal(s)
  955. char *s;
  956. {
  957.      printf("%s\n",s);
  958.      exit(-1);
  959. }
  960.